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We show that a large class of pulse coupled oscillators converge with high probability from random 
initial conditions on a large class of graphs with time delays. Our analysis combines previous local 
convergence results, probabilistic network analysis, and a new classification scheme for Type II phase 
response curves to produce rigorous lower bounds for convergence probabilities based on network 
density. These bounds are then used to develop a simple, fast and rigorous computational analytic 
technique. These results suggest new methods for the analysis of pulse coupled oscillators, and 
provide new insights into the operation of biological Type II phase response curves and also the 
design of decentralized and minimal clock synchronization schemes in sensor nets. 



Synchrony in systems of pulse-coupled oscillators 
(PCOs) is an important feature in physics, biology and 
engineering. Synchronization can range from being a 
pathological breakdown, as in epilepsy Pj to one of vi- 
tal importance, such as in the proper functioning of the 
heart's sinoatrial node [H [3], to a framework to under- 
stand complex systems [H [S] . Additionally, there are at- 
tempts to utilize the simplicity of PCO synchronization 
to synchronize wireless sensor networks [6]-E]- However, 
many of the idealized models inspired by synchroniza- 
tion are not able to synchronize when the system has a 
complicated graph structure and time delays - aspects 
expected in real physical systems. In order to deal with 
these issues, previous studies have considered oscillators 
augmented with memory [71 [10] , infinite spatial density 
[To] or indegree normalization [51 [10] . While these stud- 
ies have shown linear stability [S] , or other forms of local 
convergence [3 [10] , global convergence in these settings 
has either been shown to be impossible [^ or remains 
unknown. 

Alternatively, a class of oscillators with Type II phase 
response curves (PRCs), have been connected to syn- 
chronizing behavior theoretically [111 I12j and in nature 
[H [131 [Hj . The distinguishing feature of oscillators with 
Type II PRCs is that an oscillator's phase can either be 
decreased (inhibited) or increased (excited), depending 
on the internal state of the oscillator. In this paper we 
focus on PCOs with a particular class of type II phase 
response curves, introduced in our previous paper [T5] . 
which resembles those in nature [21[I3] and are well suited 
for handling complex topologies and time delay. We also 
show how leveraging the main theorem from [15 allows 
for a computational analytic routine yielding a fast and 
rigorous estimate of the convergence probability of a sys- 
tem of PCOs. Furthermore, we provide rigorous lower 
bounds that guarantee the performance of this compu- 
tational analytic approach and display how the proba- 
bility of synchrony converges to 1 in highly connected 
graphs. This result is of biological relevance to the sit- 
uations where synchrony is brought about via Type II 



PRCs, and is a useful guide for the construction of PCOs 
in sensor nets. 

Previous work found that a class of Type II PRC, 
denoted "Stong Type 11" or "STII" (described later), 
could consistently converge to synchrony on fairly com- 
plex graphs with time delays [15j . This convergence was 
explained by showing that these PCOs would converge to 
synchrony if their phases were inside a critical range po , 
essentially showing an l°° ball of stability. This showed 
that with well-tuned parameters the system is robust to 
any individual oscillator error or a combination of small 
errors; explaining the possibility of synchrony, but not 
the ubiquity of it in numerical simulations. For example, 
if the critical range is ^ of the phase interval, then the 
probability that a system of n oscillators with uniform 
random initial conditions starts in the critical regime is 
i , which is exponentially small in the system size; how- 
ever, numerical experiments show that convergence is in 
fact highly likely and our analysis explains this. In par- 
ticular, we use network analysis to expand on the local 
understanding of stability, showing that if node indegrees 
follow a simple scaling law, random initial conditions in 
any size system are very likely to collapse to the critical 
range. 

This analysis of STII oscillators sheds light on some of 
the natural questions regarding Type II phase response 
curves, and the contribution of excitation and inhibition 
to synchronization in Type II PRCs. For example, it's 
clear that the important aspect of the excitatory end of 
a Type II phase response curve is that it allows for fir- 
ing cascades, yet previous analytic results have tended 
to focus on the importance of inhibition when the sys- 
tem has time delays [S] [TSl [T5] . In contrast the result in 
this paper classifies the excitation in a type II PRC into 
different discrete classes, each class corresponding to its 
ability to cascade and a lower bound on the probability 
convergence. 

To understand the strong convergence of STII oscilla- 
tors, consider the following PCO model: there are n os- 
cillators on a strongly connected aperiodic directed graph 
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G. Each oscillator j's state is described by phase variable 
4>i G [0, 1], which evolves with natural frequency d(f>i/dt = 
1. When 4>i(t) = 1 the oscillator emits a pulse and its 
phase is reset from 1 to 0. This pulse is received by all 
of I's successors, S{i), time r < 0.5 later. When a pulse 
is received, oscillators process this pulse via the phase 
response curve fij{4>i)j where — ^ max(0, 0j + fij{4>j))- 
This setup for a PCO can be made to accommodate many 
of the different types of PCOs examined in the litera- 
ture. For example, the popular Mirrolo and Strogatz 
model, which uses a charging curve V , can be described 
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The central goal of this paper is to understand Pf{G), 
the proportion of phase space that converges to syn- 
chrony for a graph G and PRC /. In particular we are 
interested STII phase response curves defined similarly 
as in 1151: 



< — min(a;, t + k) : x < B 
> 0, :x>B 



where B e (0, 1) and k > are parameters. In [15] 
it was shown that if every is a STII PRC, and at 
some time all oscillator phases are within a critical range 
Po < mm{B — T, 1 — B + t), then the system will converge 
to synchrony in finite time. 

To demonstrate the ability of STII oscillators to reach 
synchrony on complex graphs with time delays, consider 
Figure [T] which shows the maximum differences between 
oscillator phases as a system is integrated for different 
PRCs and random initial conditions. Notice that not 
only are STII curves the only curves that converge, but 
for most runs, STII curves fall within the critical range 
in a single time step, despite the fact that the size of the 
critical range is exponentially small in probability space. 
As will be shown, a large portion of the basin of attrac- 
tion of synchrony in STII oscillators can be described by 
this rapid convergence to the critical range. Furthermore, 
this convergence to the critical range arises from a fun- 
damentally different mechanism, and relies on different 
properties of the PRC than the convergence inside the 
critical range. 

To understand this basin analytically, consider first, 
the most extreme STII PRC, the "strong firing" (SF) 
PRC where f^^{x) = -x ior x < B and f^^{x) = \-x 
otherwise. Notice, the response f^^ gives to any sig- 
nal causes an oscillator to have phase 0, where signals 
received after B also cause it to fire. This brief char- 
acterization of the SF PRC allows for a quick analytic 
description of one way in which the SF PRC converges 
rapidly to the critical range. 

The key insight is that if for every SF oscillator z, j re- 
ceives a signal or fires in a small window of time, denoted 
as event Ei, then every oscillator will be reset to phase 0, 
and thus all phases will be within the critical range. We 
will show that the window can be as large as 1 — s, where 
s = max(i3, 1 — B + 2t). Notice, that if every oscillator 





FIG. 1. Starting from uniform random initial conditions in a 
system with 400 nodes, trajectories (mean - solid line, middle 
50% - between dotted lines) either converge to synchrony or 
not depending on the PRC. Notice that two classes discussed 
in this paper, SF and STII oscillators converge to exact syn- 
chrony in finite time while others from ChaosOS [T7| , IEEE05 
[6], and SIAM90 [4] do not (this remains true in other mea- 
sures, not shown). Furthermore, convergence of SF and STII 
is guaranteed once the maximum difi'erence is smaller than .5 
by US]. 



receives a signal or fires at some time in [r, 1 — s + r] , 
then for all z, (f)i{l — s + t) < 1 — s < B — 2t. Since 
only oscillators with phases greater than B can generate 
signals, no new signals are being sent. Thus by time t 
later all signals will have arrived at their destinations, 
giving that for all i, s + 2t) < 1 — s + t < pQ which 

puts all oscillators in the critical range with no signals 
enroute. 

Therefore, the probability that the SF system con- 
verges can be bounded by the likelihood that every oscil- 
lator receives a signal in a window of size 1 — s time, 
P(nfi?i), which can then be bounded by using node 
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degrees, di. If at time 0, 4'j{Q) € [s, 1] then all the 
successors of j will necessarily receive a signal in time 
[r, 1 — s + r] (since s > _B, j is in the excitatory regime) 
and if (/>j (0) G [s — t,1 — t] then j must fire in [t, 1 — s + r] . 
(Note: that this is now very similar to the probability 
that a random subset of the graph will dominate it, con- 
necting it to some sensor net protocols used to find a 
Connected Dominating Set [lll[T^ and the study of dom- 
inating sets in general |20|V Thus, for uniform random 
initial conditions a simple bound on the probability of 
any oscillator i receiving a signal or firing in [t,1 — s + t], 
is simply the complement of all i's predecessors having 
phases in [0, s) and 0i £ [0, s — r] U [1 — r, 1]; yielding 
P{E^) > 1 - s'^^+\ 

Notice, that the probability of a node i failing to re- 
ceive a signal in the 1 — s time window is exponen- 
tially small in that node's indegree. These probabilities 
can be aggregated using the Union Bound, giving that: 
P(Uf£f) < Sfs'^'+i and thus PinfE,) > 1 - E^s'^'+i. 
Alternatively, a slightly stronger bound can be found us- 
ing the fact that each Ei is positively correlated, or that 
the number of nodes dominated is a submodular function 
of node subsets, giving the probability of convergence, 
Psf{G) > P{nfE,) > nf (l-s'^'+i). Thus the determin- 
istic bound from |T5] has been used to create a statement 
about convergence from random initial conditions. 

This result immediately gives a number of interest- 
ing corollaries. For example, let 6n(j}) be the mini- 
mum indegree such that the Psf{G) > p then 6n{p) < 
ln(l — p)/ ln(s) — ln(n) / ln(s), which is logarithmic in the 
system size. To give a sense of the constants, the mini- 
mum value for s occurs when B = s = .5 + t, and thus to 
ensure a 95% convergence rate in a systems with a time 
delay 5% of the period: (5„(0.95) < 5.02 -I- 1.68 ln(n); a 
result that holds for any n. 

This result can also be used to make statements about 
the convergence of SF oscillators on random graphs. Take 
an Erdos-Renyi random graph, G{n,p), where edges are 
created with independent probability p. In this case 
an application of Chernoff's inequality shows that if 
P ^ '"^"^ g(s,7) for a function g of s and some num- 
ber 7 > then as n — )■ oo, the probability of synchrony 
PsriG) ^ (1 - l/n)ei/"'"^ 1. Notice, that this re- 
quirement on p is only a constant multiple of that re- 
quired for G{n,p) to be connected, which asymptotically, 
occurs when p > {I + e)^- Thus, the degree require- 
ments grow reasonably with n, and furthermore, since 
convergence requires connectedness, our rigorous bound 
is a constant factor approximation of the actual required 
degree. 

Similarly, one can show asymptotic bounds on random 
geometric graphs, constructed by positioning nodes uni- 
formly at random on the unit d dimensional torus and 
connecting any nodes within some radius r. If r is cho- 
sen so that the expected degree r'^nO = cln(n), (where 9 
is the volume of a d dimensional unit ball) then utilizing 



results describing the minimum degree in random geo- 
metric graphs, |21! shows that the system will converge 
to synchrony as n — > oo so long as c is the greater of 
the solutions to: - = I -\ — ^ i^ln(^^). Again, 

c cms cms ^cms^ ^ ' 

a constant factor of logarithmic growth in the expected 
degree gives convergence guarantees. 

Thus far, we have shown analytically that sufficiently 
dense systems of SF PCOs will, with high probability, 
converge in a single time step to the critical range, and 
consequently, will converge to exact synchrony in finite 
time. However, synchrony is a real phenomenon, and the 
constant in our bound may be important in certain appli- 
cations. Fortunately, the analysis of the lower bound can 
be extended in a simple "computational-analytic" man- 
ner, providing a rigorous computational assisted bound. 
Figure [2] shows the relationship between indegree, the an- 
alytic lower bound, the numerical bound and two compu- 
tational analytic bounds. The first, simulates the system 
for one period and then applies the deterministic conver- 
gence result. This is fast, works surprising well and is 
far more computationally efficient than fully integrating 
the system to convergence (a time that would scale with a 
graph's aperiodic diameter). For situations where precise 
estimates of Psf{G) are important one uses the computa- 
tional analytic approach by running multiple single time 
step Monte Carlo trials and checking if the phases fall 
within the critical range. In a graph with m edges such 
a routine can be implemented by an event based simu- 
lation, and thus can run in 0(m log m) time. Whereas 
typically integration time scales with system size, our an- 
alytic bound guarantees that this computational analytic 
routine remains viable as the system size increases. 

We now consider more general PRCs, showing that 
many STII PRCs will also have Sn{p) — 0(log(n)) and 
thus will also have corresponding guarantees for random 
graph models and computational analytic routines. The 
key feature of the SF PRC was that a single signal causes 
an oscillator to reset or fire. The arguments made for an 
SF oscillator can be modified to allow for oscillators that 
require multiple signals to reset or fire. Consider the sub 
class STIIk.Tj which, as opposed to requiring 1 signal, 
will require receiving at least k signals within 1 — s — 77 
time to reset or fire. Figure [3] displays several different 
PRCs from STIIk^ for different k and 77. 

If for independent initial conditions, P(0i(O) £ [0,s + 
1])) = qi, then the probability that an oscillator i has less 
than k neighbors ready to fire is the sum of the binomial 
distribution B{di, qi) from to fc— 1. Using a well known 
bound based on Hoeffding's inequality yields that the 
probability that i has at less than k neighbors ready to 

fire is less than e ^^''i . Taking the complement 
and using the same results of independence gives that 
the probability of the system synchronizing is, 

Pf.„(G)>n^i(l-e ). 
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FIG. 2. Probability of convergence for a 400 node random 
geometric graph as a function of radius for SF (red), STIIifi 
(blue) and a STIIta (black) PCOs. Numerical results (solid) 
suggest that all three oscillators systems transition to syn- 
chrony at the same value of r. Dotted lines show an analytic 
lower bound and dashed lines show the numerical single time 
step bound. 




FIG. 3. Placing a sawtooth function underneath the exci- 
tatory portion of an STII PRC provides an upper bound on 
the number of excitations that will cause the curve to fire, 
whereas the scale of a curve's inhibition is measured in pro- 
portion to B. 



Alternatively, let c.„ = ln(n)— ln(l— p), then the system 
will converge with at least probability p if for all i, the 
expected number of firing neighbors, diqi > k — 1 + Cn + 
\/c^^ 4(fc — l)c„. Since, for fixed fc this result also scales 
0(ln(n)) then the random graph results in the SF case 
have analogs of the same order: p = 0( '"^"'' ) for Erdos- 
Renyi and r — 0{ '^"^^"j ) for random geometric graphs. 

Determining if a phase response curve is a member of 
STIIk.rj involves two steps: first, classifying the strength 
of the inhibitory section, and second, the strength of the 
excitation. 

We say that an oscillator i is ^.-inhibitory if receiving 
h signals in the inhibitory region over some span of time 
[tQ,tQ + s'], s' — 1 — s — 7/ forces (/)j(io + s') < s' . For exam- 
ple, a sufficient condition for a PRC / to be /i-inhibitory 
is if f{x) < — min(^,a;), for x < B. For such an /, h 
signals causes the oscillator either to be reset to 0, or to 
be inhibited by at least B, giving that (j){tQ + s') < s' . 

Similarly, an oscillator i is {k — h + l)-excitatory if 
receiving k — h + 1 signals in the excitatory region, in 
some time [tg , to + s'] , forces an oscillator to fire before 
tQ + s' + r]. As seen in figure [Sj a sufficient condition for 
{k — h + l)-excitability is that the PRC is greater than a 
saw tooth with slope —1 when x d [B,l — i]]. 

If an oscillator is both /i-inhibitory and {k ~ h + 1)- 
excitatory then it is a member of STIIj^ j^. These results 
can also be used to guarantee the performance of a sim- 
ilar computational analytic routine. The performance of 
a STII^fi and a STII4 Q as well as their analytic guar- 
antees can be seen in figure [2] 

Finally, it is worth noting that in many systems, such 



as systems of neurons, edges are often weighted and the 
impact that different neighbors have varies drastically 
5J. The results in this paper also extend to a weighted 
version, where each edge has weight Wij and weights are 
interpreted by the formula: fij{x) = mayi(—Wij,fij) for 
X < B and fij{x) = miii(wij, fij) for x > B, where Wij 
acts as a constraint of the phase response curve. The 
above formula for the Pf{G) remains true so long as for 
each i, Wj_i > t and for each node i there are di nodes 
j such that fij e Fh^k- 

In such a case, if k increases as 0(ln(n)) or less then 
so does 5n{p). Furthermore, if /c — !■ 00 then the require- 
ments on the phase response curve shrink to simply re- 
quiring that /'(O) = —1, and f{x) < Q ioi x < B and 
,f{x) > for X > B and that / is continuous everywhere 
except B where limj._^B- < ~£ and \\m^^B+ > e. Thus 
for very large systems our results show convergence for a 
very general class of type II oscillators. However, when 
comparing to results for "weakly coupled oscillators" one 
should recall the slightly different requirement, that in 
those cases, Wij = e and the weights are multiplica- 
tive: /ij(0) = Wijf{(j)). 

In summary, we have shown how the local convergence 
of STII pulse coupled oscillators to synchrony can be ex- 
tended probabilistically, relating graph density and phase 
response curve structure to a rigorous lower bound on the 
probability of convergence. Applying this lower bound to 
random graph models shows that the expected node de- 
gree beyond which synchronization is very likely is a con- 
stant multiple of the percolation threshold. Therefore a 
computational analytic scheme that simply sampled sin- 
gle time steps is a constant factor approximation to a 
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sampling routine that integrated for infinite time, 
extension with edge weights was also discussed. 
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